The Randomized Hough Transform Used For Ellipse Detection
A. Schuler
Abstract – The Hough transform is basically just another
integral transform such as the Fourier transform. It can be used to detect
primitive shapes such as a line in a picture. To describe a line we only need
two parameters (whether we take a and b, where a is the
slope and b the offset on the y axis, or the parameters p,
distance from origin, and Q, angle; doesn’t really matter). Now this transform
is applied (usually to a binary image which has been preprocessed with some edge
detection algorithm) to each point in the source image; the transform of each
point (x,y space) is written in the transformed space (p, Q space). By examination of the transformed space, we
are able to detect lines in the source image; because a line in the source image
corresponds to a point in the transformed space.
We can now extend this
simple approach to more complex shapes such as an ellipse. The general equation
of an ellipse contains 5 parameters, therefore our transformed space would
become 5 dimensional (since it was two dimensional for a line), what can cause a
lot of problems even for a computer to handle.
To overcome this
performance issues the RHT (Randomized Hough Transform) was introduced in 1990
by Xu et al[3]. Although the algorithm is not really applicable to curves that
involve non linear parameters, it still is possible to use the RHT for ellipse
detection if the process can be reduced to a linear problem (this has been done
by McLaughlin in 1997[2]).
In this project the author tried to implement McLaughlin’s algorithm in MATLAB. The image processing toolbox is only required for preprocessing images not the RHT itself.
Introduction
The Hough transform has been
proposed in 1962 by P.V.C Hough. It was first used for pattern recognition in
the early 1980s.
The basic concept of the
Hough Transform is quite simple, lets look at the
integral:
This equation tells us, that
each pixel (I(x,y)) is transformed to a sine wave in the parameter space. These
waves intersect with each other if the pixels that generated them lie on a line
in the source image (see picture). It should be clear now, that it is easily
possible, to extend this equation for curves of a higher order. For the example
of an ellipse, which is given by the following equation:
the
parameter space would become five dimensional which would make postprocessing of
the parameter space almost impossible (single points in the parameterspace would
now correspond to ellipses instead of lines).
McLaughlin introduced an improved version of the RHT
(Randomized Hough Transfrom) to simplify the detection
process.
Application
The question the author
asked himself was: is it possible to detect ellipses of any size or rotation in
a binary image even if they are slightly overlapping? There are many
applications where this algorithm could be used; for example the counting of
blood cells in a sample (the algorithm could also handle overlapping cells) or
even an accurate measuring of ellipses in any picture (it would be possible to
compute the exact size and location of each ellipse). Since the algorithm
returns the parameter value for each ellipse found it would be very easy to
postprocess and store (once the ellipses are found the data could be stored
additionally to the image data) the data (for instance elimination of false
detected ellipses).
Detailed description of the algorithm
The algorithm is basically a
stochastic process, what means that the outcome doesn’t have to be the same each
time. The image has to be preprocessed with an edge finding algorithm (e.g.
Sobel) and contains only binary values, i.e. pixels on an edge are 1 (=white)
and all other pixels are 0 (=black).
1. Take any 3 points out of the image that lie on an edge. The maximum and minimum distance two of the three points can have is restricted.
2. Try to find an estimate for
the center of the ellipse. This is done by the following algorithm:
Find an
estimate for the tangent through each point p1, p2 and p3 (look at neighboring
points an minimize the sum of the squared distances from these points to the
tangent) => tan1, tan2 and tan3.
Find intersection of tan1 and tan2 =>
t1, and find intersection of tan2 and tan3 => t2.
Find midpoints of p1p2
=> m1 and p2p3 => m2.
The estimated center now lies on the intersection
of t1m1 and t2m2.
If center cannot be found, go to
1
Figure 2 Finding the center of an
ellipse 1 Figure 3 Finding the center of an
ellipse 2
3.
Now the center is found
translate the points p1 .. p3 to the origin (0,0).
4.
Now compute the
parameters a, b and c since we know that 
by solving the
following equation (p1=x1,y1 p2=x2,y2
p3=x3,y3):
If there doesn’t exist a solution, go to
1
5.
Check inequality 
, which has to be true for a valid ellipse. If false, go to
1.
6.
Now compute the
parameters p, q, r1, r2 and Q since this parameter space is better behaved. Those
parameters satisfy the equation:
7.
Now search the list of
already found parameter sets (originally a list, as proposed by McLaughlin[2]
but in MATLAB one better uses a matrix). If there is a parameter set whose
values are close to the ones newly found (determined by an epsilon), then
average those two sets and increase the count of this parameter set. Each
parameter set has a so called count value, which tells how many times these same
parameters have been found.
If there is no similar set, simply add the new
parameter set to the list and set its count to 1.
8.
After a specified
amount of those parameter sets is found (so called epoch), analyze the list (or
matrix) of parameter sets. Take every set that has a count value above some
threshold and check it.
By checking we mean going back to the binary source
image and counting how many points actually lie on the ellipse specified by the
parameters found. This counting has to be done with some tolerance of course and
the amount of pixels is then normalized with r1 and r2 (to be accurate: the
amount of pixels found on the ellipse is divided by the number of pixels lying
on the rectangle defined by 2*r1 and 2*r2). This normalized value is then
compared to some threshold and the ellipse is eventually considered
existent.
9.
If there were some
ellipses found then those pixels are erased from the original image (i.e. the
image gets simplified). Now we clear the list (or matrix) of parameter sets and
start with another epoch. (i.e. go to 1)
From this description it
should be clear, that there are many parameters in this algorithm that have an
impact on the performance:
·
What’s the
maximum (minimum) distance allowed fro two points out of the three
sample points. If we set the maximum to a very high value, it will become almost
impossible to detect small ellipses (and vice versa of course). The minimum
boundary simply should the algorithm prevent from taking two points too close to
each other (so we would get a pretty bad estimate for the ellipse center and
therefore for the rest of the ellipse parameters.
Value in MATLAB rht.m (in
pixel):
minPdist
maxPdist
·
How many
parameter sets are collected in
one epoch? This value is more ore less a function of the image
size and the number of edge pixels: the bigger an image is (or the more edge
pixels are present) the more parameter samples should be taken so we get at
least one ore two sets with a count higher than 1.
Value in MATLAB
rht.m:
epoch
·
What should the
count threshold for a parameter set be? After each epoch, when we
examine the parameter set list, we need to decide which ellipse parameters are
worth checking. The higher this threshold value is, the harder will it become to
detect an ellipse. A high value also implies more set samples per epoch.
Practically a value of 2 (i.e. all parameter sets witch a count higher than 2)
has been proven quite reasonable.
Value in MATLAB rht.m:
countThreshold
·
In finding the
tangents, how big should the window around the point
examined be? To find a tangent for one of the three sample points, we need to
consider its neighboring points; the more neighboring points we consider, the
better will the estimate for the tangent be. On the other hand, the bigger this
window, the harder will it become to detect small ellipses.
Value in MATLAB
rht.m (in pixel):
tangentSearchWindowRadius
·
If we’re checking an
ellipse for its existence, we basically count all pixels lying on the ellipse
described by the parameter set examined, this cannot work without a certain
tolerance, i.e. we consider pixels lying near the
ellipse, too. This value can be set in pixels (i.e. number of pixels
that the ellipse radius can be bigger or smaller). Of course it can increase the
number of false alarms, i.e. it detects ellipses where there are none; a too
small value on the other hand makes it almost impossible to detect any ellipse
since in natural pictures the ellipses are never perfect (gaps,
irregularities…).
Value in MATLAB rht.m (in pixel):
ellipseRangeTolerance
·
The above point leads
to another issue: how to decide if the ellipse found exist. To accomplish that,
we need to normalize to pixels counted along the ellipse some how.
In this implementation it is done be dividing the number of pixels found by the
number of pixels lying on the rectangle specified by 2*r1, 2*r2. The value is
the compared to a certain threshold after which the ellipse is considered
existent or not.
Value in MATLAB rht.m (in percent, e.g. 25 means
25%):
minCoverage
·
How is a
parameter set considered ‘close’ to another set? This is very
important when we have found a new parameter set and try to add it to our list.
We have to compare it to every list member and decide, if there is already a set
with values close enough (-> threshold). If this threshold is
too small, we will always detect new parameter sets and never increase the count
of one set, this prevents detection. If this threshold is too big, we might melt
to actually different (in the source image) ellipses together (ellipses that are
spatially close).
Another consideration worth is the value for Q in our parameter set. Since were calculating the sum
of absolute errors for each parameter set in the list (i.e. errori =
abs(pi-p)+abs(q-q)+abs(r1i-r1)+abs(r2i-r2)+abs(Qi-Q), where errori is the error of parameter
set i against the set examined) and we know that the value for Q is in radians (=> small numbers!) we have to
weight this value by a factor. The other parameters don’t have to be weighted
since they are all in pixels.
Value in MATLAB rht.m:
maxDiffError
Tests
Figure 4 Test image a, 4 ellipses
no overlap, 256x256
Figure 5 Output image, input test
a
Figure 6 Test image b, 6 ellipses, overlap clipping, 256x256
For the algorithm it took
approximately equal amount of epoch to detect all ellipses. Hard to spot was
especially the very flat ellipse in the upper right corner (this was also an
issue in test image a). Less surprising is the performance for the overlapping
ellipses: almost without any false alarm they could be detected easily. The
following parameters were used:
Figure 7 Output for test image
b
We slightly change the test
image b, to test the real power of the RHT:
The image now contains
incomplete ellipses (from overlapping) as it could occur in a natural
image.
Figure 8 slightly changed test
image b, incomplete ellipses
The parameters of the
algorithm remains unchanged and so was the result. The flat ellipse still caused
the most trouble. The output could look like this (the output image is
superimposed on the input image):
Figure 9 Output to modified test
image b
c.) complex image
containing 17 ellipses, overlapping, incomplete
Figure 10 Test image c, 17
ellipses, overlapping, incomplete
This image now could be
acquired from a natural image showing some blood cells (although it’s
artificial). It contains 17 ellipses some of them overlapping. Most other
counting algorithms would now detect 14 objects (maybe with intensive
preprocessing the lower two ellipses could be separated…).
The only parameter that was
changed from test b to test c is MaxPdist (50) since we’re mostly dealing with
quite small ellipses in this test.
The algorithm didn’t show
any problems with this test set, it could always detect all ellipses (without
false alarms) in about 10 epochs.
Sample
output:
Columns 1 through 7
112.9360 84.8508 42.2839 37.0221 229.8996 202.9475 170.1096
117.8678 170.4458 223.2040 26.9583 179.6400 15.3995 146.3641
15.1776 15.0324 23.3171 18.2418 20.0866 19.3843 19.2995
21.9421 24.6337 16.5942 12.9795 17.0011 13.0494 13.2101
-0.0306 0.0346 0.6060 -0.0335 -0.1398 0.5956 -0.6915
Columns 8 through 14
100.1033 200.2431 144.6983 54.9993 31.2877 197.9204 209.5656
32.7699 74.0598 200.1267 82.8938 152.2040 223.3853 102.6654
20.2160 16.2228 20.1443 19.0936 27.6900 18.6038 22.2431
13.3878 26.3834 17.1148 12.2795 20.1083 16.7905 17.3532
0.0019 -0.4706 0.6066 0.6520 -0.6377 0.6551 0.0415
Columns 15 through 17
145.6133 117.3194 22.5593
75.0058 215.1759 119.1884
13.8327 21.3291 16.7653
11.7336 15.1020 23.5869
0.0995 -0.2305 -0.2200
Figure 11 output to test image
c
d.) natural image
containing about 10 detectable
ellipses,
Figure 12 natural image (blood
cells)
The picture has to be
preprocessed first (preProcess.m):
·
Image normalized
(values between 0 and 1)
·
Contrast
stretch
·
Noise removal by
applying a wiener filter (which gives better results than a simple median
filter)
·
Edge detection using
the Sobel operator
This results in this
image:
Figure 13 preprocessed test image
d
For detection we can use the
following parameter:
We only changed the
MinCoverage value because the cells don’t have perfect elliptic shape, so we
introduce a slightly bigger tolerance.
Again the algorithm needs
about 10 epochs to detect all ellipses. Problems can be caused by the incomplete
ellipse on the left side and the deformed ellipse in the center of the image.
Nevertheless the RHT finds almost every time all ellipses correctly. A sample
output:
Columns 1 through 7
270.3021 63.6827 336.9685 110.3143 164.8037 282.6283 177.4440
108.1948 121.0237 88.5797 34.1158 59.3212 206.4743 98.1646
24.5916 24.2187 24.8545 24.1066 21.3167 18.1552 25.1232
24.0916 21.4756 22.3462 21.0549 17.3210 20.5081 18.5974
-0.7254 0.1735 -0.5263 -0.2686 0.2958 -0.2378 0.1678
Columns 8 through 10
160.3133 349.6957 1.2822
144.1991 13.2230 96.7674
28.9921 29.1964 14.4132
26.4200 24.5893 20.5341
0.6404 -0.6144 -0.2035
And the superimposed
image:
Figure 14 output to test image d,
blood cells
Results and conclusion
The RHT performs very well
with artificial images as long as the ellipse’s radii are not too different
(results in flat ellipses). Even with the natural image the algorithm still
works pretty good. Especially the accuracy of the parameters found is quite
astonishing if one compares the input and output images. This implies that this
method of pattern recognition is suited very well for highly accurate
measurements (e.g. exact location, rotation, locus, area etc…). Since the
algorithm provides all the results in matrix form, it’s very easy to postprocess
this data.
There are some major
drawbacks though: since the algorithm works on a stochastic base, we don’t get
the same result each time we apply it to an image; even worse, it can detect
some nonexistent ellipses and/or not detect others sometimes. This behavior
suggests to run the algorithm several times over the same image and average the
results (and hopefully discard false alarms). Most images tested had a
resolution of 256x256 pixels (363x241 pixels for the natural image) and needed
about 10 epochs to find all ellipses. This process took about 1:30 Minutes on a
PIII 450MHz which is quite slow. This would make more than one run almost
impossible.
There are several reasons
for this bad performance (in terms of execution speed). The main reason is the
implementation in MATLAB, because the .m files of MATLAB are script files that
are interpreted in runtime. The author chose MATLAB because of its ability to
use easy and fast matrix calculations, what would have taken a lot of time if
implemented in c. The first version of this algorithm used a list structure as
proposed by McLaughlin, but it was so slow, that the storing of the parameters
in a matrix (although more memory consuming) was much faster than the list.
There is a lot of room for
optimization in the algorithm as implemented by the author in MATLAB. So if
somebody is interested in further pursuing this RHT he should feel free to
contact the author by email.
Bibliography
[1.] V. F. Leavers. "Shape Detection in Computer Vision Using
the Hough Transform, "
London: Springer-Verlag, (1992)
[2.] R. A.
McLaughlin, "Technical report - Randomized Hough transform: Improved ellipse
detection with comparison," Tech. Rep. 97 / 1, The University of Western
Australia, Centre for Intelligent Information Processing Systems, Dept. of E.E.
Eng., U.W.A., Stirling Hwy, Nedlands W.A. 6907, Australia, (1997).
[3.] L. Xu, E.
Oja and P. Kultanen, "A new curve detection method: Randomized hough transform
(RHT)," Pattern Recognition Letters, 11, 331-338 (1990).
Appendix, MATLAB code
Example of how to use the
code (it is assumed that the imaging toolbox for MATLAB is installed, else all
im*** functions won’t work as well as the preProcess
script):
blood1=preProcess(blood);
imshow(blood1);
[res
img]=rht(blood1,2);
imshow(
(double(img)+double(blood1))./2)
File:
rht.m
function
[res,
resImg]=rht(im, exitThresh)
%
RHT (Randomized Hough Transform)
%
This function analyzes a binary picture (pre processed with an edge
finding
%
algorithm) and tries to find ellipses of any shape and
rotation.
%
It outputs a matrix 'res' containing the 5 ellipse parameters in each
column
%
(they belong to the ellipse equation) :
%
( (y-q)*sin(theta) + (x-p)*cos(theta) )^2 / r1^2 +
%
( (y-q)*cos(theta) - (x-p)*sin(theta) )^2 / r2^2 = 1
%
where ...
%
p, q : the coordinates (x,y) of the
ellipse translated from (0,0)
%
r1 : radius of the
ellipse along the x-axis
%
r2 : radius of the
ellipse along the y axis
%
theta : angle to x axis (in radians) of ellipse rotated counter
clockwise
%
%
input : im
binary image (containing only 0 or 1)
%
exitThreh number of
unsuccessful consecutive epochs after which the
%
function terminates (this guarantees abortion of the
program)
%
output : res
matrix containing the ellipse parameters (each column
stands
%
for one single ellipse
%
resImg
image containing the ellipses found (may be used as
mask)
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
%
global variables can be accessed by every subroutine
%
contains a 6xn matrix (5 ellipse parameters and the count, i.e. how many
times
%
there has been found an ellipse in the picture with the same param
set).
global
parMat;
%
nx2 matrix containing all the coordinates (x,y) of the image that are non 0
(i.e.
%
all points that lie on an edge)
global
points;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
PARAMETER SECTION
%
%
The follwing variables may have an impact on the performance of
this
%
algorithm. They should be changed for the needs of the input image,
i.e.
%
an image containing only very small ellipses should have small
'minPdist'
%
and 'maxPdist' values (among other changes...).
%
Some variables wre global becaus they are used in sub
routines.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Value in percent.
%
This value is used by the 'checkEllipsesFound' function. The points found
%
on an ellipse are normalized with the points lieing of the rectangle
defined
%
by 2*r1, 2*r2 (i.e. points on
ellipse / points on rectangle). If this value
%
(percentage) is bigger than 'minCoverage' then the ellipse is considered
%
existent.
global
minCoverage;
minCoverage
= 30;
%
This value is used by the 'findEllipseCenter' function. It specifies a
radius
%
in pixels around a certain point. Originating in this point the algorithm
then
%
tries to match a tangent that best fits the points inside the cricle
(specified
%
by the radius).
%
A value too big may prevent detection of small ellipses, although with a
bigger
%
window size the estimate for the tanget gets better.
global
tangentSearchWindowRadius;
tangentSearchWindowRadius
= 7;
%
When adding a newly found parameter set then the difference is computed to
each
%
param set in 'parMat' (=> 5 errors each set). The sum is then taken of the
%
5 absolute values for each set. If the summed error lies below 'maxDiffError'
then
%
the two parameter sets are considered equal.
global
maxDiffError;
maxDiffError
= 4;
%
minimum and maximum distance that two of the three picked
points
%
are allowd to have (it doesn't make sense to pick points too far
away
%
from each other).
%
this parameters can affect the algorithm's ability to detect very large
%
(or very small) ellipses in the image
minPdist
= 3;
maxPdist
= 50;
%%%%%%%%%%
Quantisation is NOT used anymore in this version !!!!!!
%
Quantisation steps for each of the 5 ellipse parameters. The smaller the
values,
%
the finer will the resolution be but the harder will it be to find
more
%
than one point set (3) with equal parameters.
pQuant
= 2;
qQuant
= 2;
r1Quant
= 2;
r2Quant
= 2;
thetaQuant
= 0.2;
%
Number of loops done before parameter data is analyzed
%
i.e. algorithm picks epoch * 3 points.
%
The bigger this number, the more likely it is to find more than
one
%
point set with equal ellipse params.
epoch
= 80;
%
To consider an ellipse in the image as existent, the algorithm needs
to
%
find more than 'countThreshold' occurences of an ellipse param
set.
countThreshold
= 2;
%
Value in pixels.
%
If the algorithm checks the existence of an ellipse it looks for
the
%
points that satisfy the ellipse equation but also within a radius by
%
'ellipseRangeTolerance' pixels bigger and smaller, since usually the
points
%
of an ellipse in the binary image don't lie exactly on the locus of
the
%
ellipse described by the 5 parameters.
%
The ellipses are also drawn with this tolerance in the output
image.
ellipseRangeTolerance
= 4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
END OF PARAMETER SECTION
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
initialization of variables
points
= [];
parMat
= [];
res
= [];
%
get number of rows and columns of input image
[M,N]=size(im);
%
result image (has same size as input image)
resImg=zeros(M,N);
%
find all non zero elements (indices are stored in di,dj)
[di,
dj]=find(im);
%
points now contains all x coordinates in the first column
%
and all y coordinates in the second column
points=[dj
di];
exitCond
= 0;
%
program exits after 'exitThresh' unsuccessful epochs
while
exitCond
<= exitThresh
%
get number of points found
num=size(points,1);
%
if number is smaller than 3 -> exit
if
num
<= 3
break;
end
for
m = 1:epoch
picks =
0;
while picks < 3*num
% get three points (indices to rows of
'points')
pickPoints=round(rand(1,3).*num+0.5);
p1=points(pickPoints(1),:);
p2=points(pickPoints(2),:);
p3=points(pickPoints(3),:);
picks = picks + 1;
% analyze points if the lie within the
specified range
if checkPoints(p1, p2, p3, minPdist,
maxPdist) == 1
% estimate ellipse
center
center=findEllipseCenter(p1,
p2, p3);
if center(1) >= 0 & center(2)
>= 0
% get other three
parameters
params=findEllipseParams(p1,
p2, p3, center);
if params(1) ~= -1
% break only if we've found a valid
param set
break;
end
end
end
end
%
there obviously was no ellipse so exit the program
if picks >= 3*num
break;
end
% add the parameters found to our 'parMat'
matrix
%addParam(params,[pQuant qQuant r1Quant r2Quant
thetaQuant]);
addParam(params);
end
%
if there were no valid paramsets found -> exit program
if
isempty(parMat) == 1
break;
end
%
find parameter sets with a count higher than a certain
threshold
col=findParamAboveThresh(countThreshold);
%
extract maximum value
%[value
index]=max(parMat(6,col));
%col = col(index);
%
check those ellipses found
if
isempty(col) == 0
found =
checkEllipsesFound(col, ellipseRangeTolerance);
else
found =
[];
end
if
length(found) == 0
exitCond =
exitCond + 1;
else
%
if there were ellipses successfully detected then add those
% parameters to
the output matrix 'res'
exitCond =
0;
for n = 1:length(found)
res = [res parMat(1:5,found(n))];
end
end
disp(sprintf('End of epoch, %d ellipses
found\n',length(found)));
%
reset 'parMat' for next epoch
parMat=[];
end
%
draw the ellipses found in the output image 'resImg'
for
n
= 1:size(res,2)
resImg=drawEllipse(res(1:5,n),
resImg, ellipseRangeTolerance);
end
File
addParam.m:
function
addParam(paramSet)
%
Add a parameter set of an ellipse to the 'parMat' matrix.
%
%
input : paramSet (1x5) containing p,q,r1,r2,theta
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
parMat
%
assure that all values are real
paramSet
= real(paramSet);
%
round the parameters accoding to theri param steps
%adjParamSet=paramSet+(quant./2);
%modSet=mod(adjParamSet,
quant);
%paramSet
= adjParamSet-modSet;
%
catch empty matrix case
if
isempty(parMat)
parMat = [paramSet
1]';
else
%
check if parameter set is already present in 'parMat'
col
= findParam(paramSet);
%
if there is no col with the same parameters we simply add the
% new parameter
vector
if
col
== 0
parMat = [parMat
[paramSet 1]'];
count =
1;
else
%
there already is a same parameter vector, so the count is increased by
1
% if there already is a similar parameter set then
the old set (weighted by
% its count) is averaged with the new param set
(weight 1)
count =
parMat(6,col);
parMat(1:5,col)
=( count*parMat(1:5,col) + paramSet' )./ (count + 1);
parMat(6,col) =
count + 1;
end
%disp(sprintf('Found valid ellipse param set: %f %f
%f %f %f %d\n',paramSet(1), paramSet(2), paramSet(3), paramSet(4), paramSet(5),
count));
end
File
checkEllipsesFound.m:
function
found=checkEllipsesFound(col,
tolerance)
%
Returns a vector containing the column indices (of
'parMat')
%
that were positively identified as ellipses.
%
%
input : col
vector containing column indices to 'parMat'
%
that eventually may describe a valid ellipse
%
tolerance tolerance in
percent (e.g. 15). Points within a radius
%
'tolerance'
bigger and smaller than the actual ellipse
%
are considered for evaluating the ellipse
%
%
output : found vector
containing column indices to 'parMat'
%
that describe a valid ellipse
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
parMat;
global
points;
found
= [];
global
minCoverage;
%
if more than the percentage epsilon points are found near the ellipse then the
ellipse is
%
considered being existent
epsilon
= minCoverage / 100.0;
uTol
= 1+tolerance/100.0;
lTol
= 1-tolerance/100.0;
for
n
= 1:length(col)
p =
parMat(1,col(n));
q =
parMat(2,col(n));
r1 =
parMat(3,col(n));
r2 =
parMat(4,col(n));
theta =
parMat(5,col(n));
tol = tolerance /
max(r1,r2);
if
tol >= 1;
tol =
0;
end
uTol = 1 +
tol;
lTol = 1 -
tol;
if
r1 > 0 & r2 > 0
el = (
(points(:,2)-q).*sin(theta) + (points(:,1)-p).*cos(theta)).^2/(r1^2)+(
(points(:,2)-q).*cos(theta)-(points(:,1)-p).*sin(theta)
).^2/(r2^2);
% get indices to points that lie on the ellipse
examined
ellipsePoints =
find(el <= uTol^2 & el >= lTol^2);
if length(ellipsePoints)/(4*(r1+r2)) >=
epsilon
found = [found col(n)];
removePoints(ellipsePoints);
end
end
end
function
removePoints(indicesVector)
%
Function removes points from 'points' matrix row indexed
by
%
'indicesVector'
%
%
input : indicesVector contains row indices to the
'points' mtarix
%
marking all points to be removed (e.g. [3 6 50 53]
%
removes rows 3,6,50 and 53 from 'points' matrix)
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
points;
newVec
= [];
vec
= [1:length(points)];
num
= length(indicesVector);
pLength
= length(points);
ind
= 1;
for
n
= 1:pLength
if
vec(n) ~= indicesVector(ind)
newVec = [newVec
vec(n)];
else
if ind
== num
break;
end
ind = ind +
1;
end
end
%
append the remaining indices
if
n
< pLength
newVec = [newVec
[n+1:pLength]];
end
%
create new points array (without the points to beremoved)
points
= points(newVec,:);
File
checkPoints.m:
function
ok=checkPoints(p1,
p2, p3, mini, maxi)
%
Check distances of points to each other. The longest
%
and the shortest side of the triangle (p1,p2,p3) are
%
examined.
%
%
input : p1 coordinates of
point one (vector x,y)
%
p2
coordinates of point two (vector x,y)
%
p3
coordinates of point three (vector x,y)
%
mini minimum
side length allowed for triangle (p1,p2,p3)
%
maxi maximum
side length allowed for triangle (p1,p2,p3)
%
%
output : ok has value 1 if
points satisfy criterias or 0 if they don't
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
%
calculate the 3 sides
a=p2-p3;
b=p1-p3;
c=p1-p2;
%
calculate norm of each side (=length)
norms=[norm(a,2)
norm(b,2) norm(c,2)];
%
discard if one side is too small or too big
if
max(norms)
> maxi
ok=0;
elseif
min(norms)
< mini
ok=0;
else
ok=1;
end
File
deleteParam.m:
function
deleteParam(paramSet)
%
Delete a parameter set of an ellipse from the 'parMat'
matrix.
%
%
input : paramSet (1x5) containing p,q,r1,r2,theta
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
parMat;
numCols
= size(parMat,2);
col
= findParam(paramSet);
if col
> 0
if
col == 1
parMat =
parMat(:,2:numCols);
elseif
col == numCols
parMat =
parMat(:,1:numCols-1);
else
parMat
= [parMat(:,1:col-1) parMat(:,col+1:numCols)];
end
end
File
drawEllipse.m:
function
outImg=drawEllipse(paramSet,
image, tolerance)
%
Draws an ellipse into the image 'image', i.e. it sets all
points
%
to 1 that satisfy the ellipse euqation given by
'paramSet'.
%
%
input : paramSet vector containing the 5
ellipse parameters (p,q,r1,r2,theta)
%
image
matrix specifying the image where the ellipse should
be
%
drawn on
%
tolerance tolerance in
percent (e.g. 15), i.e. all points that lie within
%
a radius 'tolerance' percent bigger and smaller
than
%
given by the ellipse equation ('paramSet') are set to
1
%
( = white)
%
%
output : outImg
Image containing the input image 'image' plus the
%
ellipse drawn.
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
outImg=image;
%uTol
= 1+tolerance/100.0;
%lTol
= 1-tolerance/100.0;
xDim
= size(image,2);
yDim
= size(image,1);
p
= paramSet(1);
q
= paramSet(2);
r1
= paramSet(3);
r2
= paramSet(4);
theta
= paramSet(5);
tol
= tolerance / max(r1,r2);
if
tol
>= 1;
tol = 0;
end
uTol
= 1 + tol;
lTol
= 1 - tol;
coords=[round([0:(xDim*yDim-1)]'./yDim+0.5)-1
mod([0:(xDim*yDim-1)]',yDim)];
el
= ( (coords(:,2)-q).*sin(theta) + (coords(:,1)-p).*cos(theta)).^2/(r1^2)+(
(coords(:,2)-q).*cos(theta)-(coords(:,1)-p).*sin(theta)
).^2/(r2^2);
pts
= find(el <= uTol^2 & el >= lTol^2);
for
n
= 1:size(pts,1)
outImg(coords(pts(n),2)+1,coords(pts(n),1)+1) = 1;
File
findEllipseCenter.m:
function
center=findEllipseCenter(p1,
p2, p3)
%
Estimate the center of an ellipse given three points. For
each
%
point the tanget of the ellipse is estimated (giving tang1,tang2,tang3
for
%
p1,p2,p3). The intersections of the tangents tang1,tang2 => t1 and
%
tang2,tang3 => t2 is computed. The center no lies on the intersetcion
of
%
the line originating in t1 going through the center of the side c
(p1..p2)
%
and the line originating in t2 going through the center of the side a
(p2..p3).
%
%
input : p1 coordinates of
point one (vector x,y)
%
p2
coordinates of point two (vector x,y)
%
p3
coordinates of point three (vector x,y)
%
%
output : center coordinates of the ellipse center
(vector x,y)
%
value may be [-1 -1] if not found
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
points;
global
tangentSearchWindowRadius;
%
calculate the 3 sides
a=p2-p3;
b=p1-p3;
c=p1-p2;
%
midpoint on side c
mc=p2+c.*0.5;
%
midpoint on side a
ma=p3+a.*0.5;
%
get the three tangents
tang1
= findOptimalTangent(p1, tangentSearchWindowRadius);
tang2
= findOptimalTangent(p2, tangentSearchWindowRadius);
tang3
= findOptimalTangent(p3, tangentSearchWindowRadius);
%
cacluate points where tangets intersect
%
tang1 and tang2
equ=[tang1'
-tang2' -(p1-p2)'];
equ=rref(equ);
%
tangets may be parallel at this point (eventhough very
unlikely...)
%
if matrix doesn't look like: [1 0 a]
%
[0 1 b] then tangets must be parallel
if
isequal(equ(1:2,1:2),eye(2,2))
~= 1
center = [-1
-1];
return;
end
%
intersection of tangent1 and tangent2
t1=
p1+tang1.*equ(1,3);
%
tang2 and tang3
equ=[tang2'
-tang3' -(p2-p3)'];
equ=rref(equ);
if
isequal(equ(1:2,1:2),eye(2,2))
~= 1
center = [-1
-1];
return;
end
%
intersection of tangent2 and tangent3
t2=
p2+tang2.*equ(1,3);
%
find intersection of lines t1mc and t2ma
v1=mc-t1;
v2=ma-t2;
equ=[v1'
-v2' -(t1-t2)'];
equ=rref(equ);
if
isequal(equ(1:2,1:2),eye(2,2))
~= 1
center = [-1
-1];
return;
end
center=t1+v1.*equ(1,3);
File
findEllipseParams.m:
function
params=findEllipseParams(p1,
p2, p3, o)
%
This function finds the remaining three ellipse parameters
%
if the center and three points lieing on the ellipse are
known.
%
%
input : p1 coordinates of
point one (vector x,y)
%
p2
coordinates of point two (vector x,y)
%
p3
coordinates of point three (vector x,y)
%
o
coordinates of center (vector x,y)
%
%
output : params vector containing the 5 ellipse
parameters
%
(p,q,r1,r2,theta) or (-1,-1,-1,-1,-1) if the
%
parameters don't exist for the given combination.
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
%
translate the ellipse so the center lies in [0 0]
p1t=p1-o;
p2t=p2-o;
p3t=p3-o;
%
get param p,q,a,b,c (they satisfy the ellipse equation:
%
a*(x-p)^2 + 2*b*(x-p)*(y-q) + c*(y-q)^2 = 1 )
equ=[conv(p1t,p1t)
1;conv(p2t,p2t) 1;conv(p3t,p3t) 1];
equ=rref(equ);
if
isequal(equ(1:3,1:3),
eye(3,3)) == 0
params = [-1 -1 -1 -1
-1];
else
a=equ(1,4);
b=equ(2,4);
c=equ(3,4);
%
check inequality a*c-b^2 > 0
if
(a*c)-b^2
<= 0
params = [-1 -1
-1 -1 -1];
else
if a
> 0
% write back the transformed parameters
(p,q,r1,r2,theta)
params=[o(1)
o(2) transformEllipseParams(a, b, c)];
else
params
= [-1 -1 -1 -1 -1];
end
end
end
function
tangent=findOptimalTangent(p,
radius)
%
Find the tangent that best fits the point around 'p', these points lie
within
%
a radius 'radius' from the point 'p'.
%
The fucntion uses the least squares method to find the optimal tangent, i.e.
it
%
minimizes the sum of distances squared of the points to the tangent. Usually
the
%
y value of the tangent's slope is always 1 but there are special cases that have
to
%
be considered (such as: all points are symmetric around the y axis => tangent
is
%
(1,0) ).
%
%
input : p
coordinates of the originating point for the tangent (vector
x,y)
%
radius radius in pixels that
is searched around 'p' for neighbouring points
%
%
output : tanget slope of the optimal tangent found (x,y)
originating in 'p'
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
points;
%
find points that lie within the specified radius from the point
p
radii=(points(:,1)-p(1)).^2+(points(:,2)-p(2)).^2;
inrange=find(radii
<= radius^2);
%
center all points found at [0 0]
checkpoints=[points(inrange,1)-p(1)
points(inrange,2)-p(2)];
%
suare sum of ys
ysqu=sum(checkpoints(:,2).^2);
%
square sum of xs
xsqu=sum(checkpoints(:,1).^2);
%
porduct sum
prod=sum(checkpoints(:,1).*checkpoints(:,2));
%
if all points are symmetric to either the x or y axis then
%
thier product sum can be 0 (the only solutions for a tangent can then be [1 0]
or [0 1]
if
prod
~= 0
x1=0.5*(
(-ysqu+xsqu+sqrt(ysqu^2-2*xsqu*ysqu+xsqu^2+4*prod^2))/prod);
x2=0.5*(
(-ysqu+xsqu-sqrt(ysqu^2-2*xsqu*ysqu+xsqu^2+4*prod^2))/prod);
y1=1;
y2=1;
else
x1=1;
x2=0;
y1=0;
y2=1;
end
%
[x1 y1] and [x2 y2] are now perpendicular, to find the solution we have to
calculate
%
the sum of distances for both solutions
d1=getDistSum([x1
y1], checkpoints);
d2=getDistSum([x2
y2], checkpoints);
if
d1
< d2
tangent=[x1
y1];
else
tangent=[x2
y2];
end
function
col=findParam(paramSet)
%
Find a specified parameter set of an ellipse in the 'parMat'
matrix.
%
%
input : paramSet (1x5) containing
p,q,r1,r2,theta
%
%
output : col
column index of the set found in 'parMat' or 0 if the
set
%
couldn't be found.
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
parMat;
global
maxDiffError;
thetaWeight
= 5;
parSize=size(parMat);
searchMat=[paramSet(1)*ones(1,parSize(2));
paramSet(2)*ones(1,parSize(2)); paramSet(3)*ones(1,parSize(2));
paramSet(4)*ones(1,parSize(2));
(thetaWeight*paramSet(5))*ones(1,parSize(2))];
parMat2
= [parMat(1:4,:);parMat(5,:).*thetaWeight];
diffMat=parMat2
- searchMat;
sumVec=sum(abs(diffMat),1);
%col
= find(sumVec == 0);
%
find all paramsets where their arror sum is below a certain
threshold
col
= find(sumVec < maxDiffError);
if
isempty(col)
== 0
%
if there are paramsets with a low error sum then take out the one
with
%
the smallest error sum
[value col] =
min(sumVec);
end
if
isempty(col)
col = 0;
end
function
col=findParamAboveThresh(thresh)
%
Find the parameter sets in the 'parMat' matrix that have a count
%
higher than 'thresh'.
%
%
input : thresh value specifying the threshold for
detection, i.e. the count
%
of a param set must be bigger than thresh to be
detected.
%
%
output : col
column indices of the sets found in 'parMat'
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
global
parMat;
function
distSum=getDistSum(tangent,
pointSet)
%
Caclulates the distances of each point to the tangent
(located
%
in (0,0) ) and return its sum
%
%
input : tanget vetor specifying the
tangent's slope (vector x,y)
%
pointsSet 2xn matrix containing the points to be examined
(x,y)
%
%
output : distSum sum of distances of the points
given to the tangent
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
%
get number of points to examine
numPoints
= size(pointSet,1);
tmat=ones(numPoints,3);
tmat=[tmat(:,1).*tangent(1)
tmat(:,2).*tangent(2) tmat(:,3).*0];
%
append a third column to the 'pointSet' because the crocc product can
%
only be taken from 3d vectors
pmat=[pointSet
zeros(numPoints,1)];
if
numPoints
<= 1
crossP =
cross(pmat,tmat);
else
crossP
= cross(pmat,tmat,2);
end
%
dist contains distances from every point to the line
defined
%
by the tangent vector
%
avoid zero division
if
norm(tangent,2)
> 0
dist=(crossP(:,3)./norm(tangent,2));
%
take sum of absolute values
distSum=sum(abs(dist));
else
%
if tangent vecor was (0,0) return arbitrary high value
distSum
= 100;
end
function
numPixels=getEllipse(paramSet,
tolerance)
%
Returns the number of pixels satisfying the ellipse
%
equation fiven by 'paramSet'.
%
Points count towards 'numPixels' value if they lie within
%
a radius 'tolerance' percent bigger and smaller than the
ellipse
%
described by 'paramSet'.
%
%
input : paramSet vector containing the 5
ellipse parameters (p,q,r1,r2,theta)
%
tolerance tolerance in
percent (e.g. 15)
%
%
output : numPixels number of
pixels lieing on the ellipse
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
if
(paramSet(3)
~= 0) & (paramSet(4) ~= 0)
xDim =
round(paramSet(3))+1;
yDim =
round(paramSet(4))+1;
coords=[round([0:(xDim*yDim-1)]'./yDim+0.5)-1
mod([0:(xDim*yDim-1)]',yDim)];
ellipse =
(coords(:,1).^2)./(paramSet(3)^2) +
(coords(:,2).^2)./(paramSet(4)^2);
%
find points within ellipse
pts
= find(ellipse <= (1+tolerance/100) & ellipse >=
(1-tolerance/100));
%
number of points covered by ellipse (multiplied by 4 since we
took
% the points only in the first
quadrant
numPixels
= 4*length(pts)-4;
else
numPixels
= 0;
end
function
params=transformEllipseParams(a,
b, c)
%
Transform the parameters of the ellipse equation:
%
a*x^2 + 2*b*x*y + c*y^2 = 1
%
to parameters of the ellipse equation:
%
(y*sin(theta)+x*cos(theta))^2/r1^2 + (y*cos(theta)-x*sin(theta))^2/r2^2 =
1
%
%
input : a value for
parameter a
%
b
value for parameter b
%
c
value for parameter c
%
%
output : params vector (3x1)
containing the transformed parameters r1,r2,theta
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
if
a
< 0 & c < 0
a = abs(a);
c = abs(c);
end
r2=1/2*sqrt(2)*sqrt(-(-b^2+c*a)*(-c-a+sqrt(c^2-2*c*a+a^2+4*b^2)))/(-b^2+c*a);
r1=1/2*sqrt(2)*sqrt((-b^2+c*a)*(c+a+sqrt(c^2-2*c*a+a^2+4*b^2)))/(-b^2+c*a);
if
b
~= 0
ok = 1;
else
ok
= 0;
end
if
r1
== r2
ok = 0;
end
if
r1
== 0 | r2 == 0
ok = 0;
end
if
ok
== 1
fac=b*(r1^2*r2^2)/(r2^2-r1^2);
y=-1/4*sqrt(2+2*sqrt(1-4*fac^2))*(-1+sqrt(1-4*fac^2))/fac;
x=1/2*sqrt(2+2*sqrt(1-4*fac^2));
theta=-i*log((x+i*y)/sqrt(x^2+y^2));
else
theta
= 0;
end
res
= [a b c] - transformEllipseParamsInv(r1, r2, theta);
if
sum(abs(res))
< 1e-4
params=[r1 r2
theta];
else
params=[r2
r1 -theta];
end
function
params=transformEllipseParamsInv(r1,
r2, theta)
%
Transform the parameters of the ellipse equation:
%
(y*sin(theta)+x*cos(theta))^2/r1^2 + (y*cos(theta)-x*sin(theta))^2/r2^2 =
1
%
to parameters of the ellipse equation:
%
a*x^2 + 2*b*x*y + c*y^2 = 1
%
%
input : r1 value for
parameter r1
%
r2
value for parameter r2
%
theta value for
parameter theta
%
%
output : params vector (3x1) containing the transformed
parameters a,b,c
%
%
by Andrew Schuler (a.schuler@usa.net) March, 2000
if
r1
> 0 & r2 > 0
a=(r1^2*sin(theta)^2+r2^2*cos(theta)^2)/(r1^2*r2^2);
b=(sin(theta)*cos(theta)*(r2^2-r1^2))/(r1^2*r2^2);
c=(r1^2*cos(theta)^2+r2^2*sin(theta)^2)/(r1^2*r2^2);
else
a=-1;
b=-1;
c=-1;
end
function
outImg=preProcess(inImage)
%
normalize the image
inImage
= im2double(inImage)./255;
%
contrast stretch the image
inImage
= cst(inImage);
%
remove noise from image osing a wiener filer
inImage
= wiener2(inImage,[10 10]);
%
enhance edges
inImage=edge(inImage,
'sobel');
outImg
= inImage;
function
out=cst(in)
%
contrast stretch
mini=min(min(in,[],1));
maxi=max(max(in,[],1));